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Using the Algebraic Bethe Ansatz (ABA) in conjunction with a simple Monte Carlo (MC) sam- 
pling technique, we study the problem of the decoherence of a central spin coupled to a nuclear 
spin bath. We describe in detail the full crossover from strong to weak external magnetic field 
field, a limit where a large non-decaying coherence factor is found. This feature is explained by 
Bose-Einstein-condensate (BEC)-like physics which also allows us to argue that the corresponding 
zero frequency peak would not be broadened by statistical or ensemble averaging. 




Using the spin of a single electron (or hole) trapped 
in semiconductor-based quantum dots has been a long 
standing proposal for a possible implementation of a 
l| for which single-spin readout and coherent con- 
is now possible. In these setups, the isotropic 
Fermi contact hyperfine interaction of the trapped cen- 
tral spin with the bath of nuclear spins present in the sub- 
strate is know to be the essential source of decoherence 
leading, over time, to a loss of the information encoded 
in any prepared state of the central spin. Denoting by 
the central spin-i, by Ij the N nuclear spins, and by g 
and gn the respective couplings to the external magnetic 
field h, the Hamiltonian reads 
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H = ghS^ +gnhY, I 3 + ^" • ^ 
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Dropping an irrelevant constant the first two terms can 
be simplified 0,01 to BS^ with B ^ [g - g^)h. 

In the regime of strong magnetic field B > A, with A 
being N times the largest A^ , the time evolution of the 
central spin {So{t)) can be described perturbatively 0,0- 
lij in the flip-flop terms IJ + h.c). In a de- 

tailed study to fourth order Coish et al. Q found that 
the Larmor precession undergoes a typical exponential 
decay which is supplemented by additional low-frequency 
envelope modulations. In the opposite limit of weak (or 
intermediate) fields a variety of methods including semi- 
classics, exact diagonalization, or the ABA have been 
applied 0, For unpolarized spin baths at i? = 0, 

extrapolating a semi-classical approach 0] to the con- 
tinuum limit, as well as time-dependent mean field the- 
ory [ll[ for very large systems with N = 16000 showed 
logarithmic decay {S§{t)) ^ l/ln(t) at long times. How- 
ever, all approaches in the weak field limit are either 
based on mean- field/semi-classics or are restricted to ei- 
ther very small system sizes < 20, specific bath polar- 
izations, or the short-time behavior. 

In this work we aim to overcome these limitations by 
developing a hybrid method based on the ABA in com- 
bination with a direct MC sampling of the exact eigen- 
states of the system, which we apply to the transverse 



relaxation of the central spin. The method provides a 
non-perturbative full quantum treatment giving us ac- 
cess to the complete crossover from strong to weak field 
without any restriction on the initial polarization of the 
bath. Due to the trivial time evolution in the eigenbasis, 
it also provides access to real-time dynamics at arbitrarily 
long times without accumulated errors. While it is only 
usable for modest system sizes, recent developments have 
drastically sped up the computation time. In this work 
we treat systems containing up to A^ = 48 nuclear spins 
whose Hilbert space is ~ 2 x 10^ times larger than the 
ones previously treated with exact methods 0, [l^l ■ 

For the dynamics of {S~^{t)) we find the following con- 
clusions: (i) At strong to intermediate fields we observe 
exponential decay with weak modulations and thus con- 
firm the perturbative picture 0- (ii) For smaller fields, 
these weak modulation become dominant leading to long- 
lived slow oscillations, (iii) As i? ^ 0, an initial rapid de- 
cay is followed by the formation of a non-decaying coher- 
ent fraction whose amplitude and phase is in one-to-one 
correspondence with the initial state of the central spin. 
This lack of decay is explained as a consequence of BEC- 
like physics, (iv) We argue that the non-decaying fraction 
survives ensemble averaging and is at most logarithmi- 
cally dependent on finite size effects. Consequently, we 
expect it to be experimentally observable. 

To be specific, we will consider nuclear spins ^ al- 
though this restriction is not imposed by integrability 
or the numerical method. The couplings are chosen as 
Aj = ^e^^-'^^^/^ corresponding to a Gaussian electronic 
wavefunction in a 2D quantum dot 0. We study in detail 
the time-evolution of an initial coherent superposition of 
the central spin: (j-f)-) + |J|)) related to the T2 trans- 
verse relaxation time. We assume the initial condition to 
be a product state I^-q) = ^ + |J|)) ® l^'bath), with 
^I^bath) describing the nuclear spins. Narrowing tech- 
niques [l^l will typically allow one to create a superposi- 
tion of (nearly-)degenerate eigenstates of the Overhauser 
operator ft,^ — J^f^i ^j^j ■ Provided there is no special 
phase relation between them, the density matrix can be 
reduced to its diagonal terms, leading to an averaged sum 
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over the represented eigenstates. Here we avoid the 
average by using a single specific zero-polarization " max- 
imal entropy" state: l^'bath) = lii t2 is t4 • • • ) with 
the eigenvalue hf^-^ = X]jLi(~l)"'^j/2- Having the up- 
pointing spins uniformly spread out over the full range of 
available coupling strengths should make this state suffi- 
ciently generic to capture the dominant features of a di- 
agonal ensemble average. We will confirm the validity of 
this particular choice by explicitly comparing randomly 
generated initial nuclear configurations. 

Our main quantity of interest is the central spin coher- 
ence factor which, projecting on the eigenbasis of iJ, can 
be written as 

= (*(<)! S+ !*(<)) = ^C„,„e'("'— )*, (2) 



where Cm,n = (it; *bath \lpm) {lpm\ |-0n) {ipn |JJ-; *bath) 

with liprn) and \tpn) denoting eigenstates of the Hamil- 
tonian ([1]) with energies LOm and uj„ respectively. Since 
Sq + J2j=i is conserved, \ipm) must contain one more 
up-spin than |V'n)- 

Any eigenstate of the central spin model ([T]) is entirely 
defined by M complex rapidities {Ai...Aj\f} which are a 
solution of the system of M coupled non-linear algebraic 
Bethe equations 
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Here ek = ~1/Ak and eo — 0. For every solution of 
([5]) the corresponding (unnormalized) eigenstate is ob- 
tained by the repeated action, for each rapidity, of a 

generalized creation operator S^(Ai) = X!fcLo x^-^ ' i-^- 
|{Ai...Am}) = nfiiS+(AO 1^;;; ... ;). The correspond- 
ing eigenenergy is then given by 



E{{Xi...Xm}) 
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Due to the level of difhculty of systematically finding 
solutions to Eq. ([3]) we solve instead for a different set 
of variables A(e,) = y^*^, — which can be shown to 
obey a simpler set of quadratic equations jU, [15[ . Any 
given solution to these equations is found starting from 
the trivial B = oo solutions where an ensemble of M 
spins are pointing up and the remaining N — M are point- 
ing down. These configurations are deformed by a step- 
wise ramping of the l/i3-parameter until the desired B 
value is reached 15| . Finally, it was recently shown that 
the scalar products matrix elements defining Cm,n in Eq. 
([2]) can be written, in terms of A(ei), as determinants of 
A^-|-l by A^-|-l matrices The resulting fast algorithm 
allows us to define a probability Pm,n = \Cm,n\ for any 
pair of eigenstates (m, n) and perform the double sum 
in Eq. ^ using a simple Metropolis algorithm (see [l^ 



for another example of combining MC with ABA) . Start- 
ing from a randomly selected pair of i? = oo eigenstates, 
we first deform them to the desired finite- i? eigenstates 
I'fpm) , \fpn) and compute the probability Pm,n, frequen- 
cies ujrn, and the sign s^^fi — sgn^C^^^). A new pair 
is then generated by randomly selecting one of the two 
B = oo configurations and minimally changing it by ex- 
changing a randomly selected pair of up and down spin. 
Deforming this new configuration to finite B we com- 
pute Prn,n' (assuming state n was modified) and accept 

the new pair (to, n') with probability max ^1, ^ . Re- 
peating the procedure generates a list of configurations 
{rria^na) distributed according to Pm^,na, such that 



< ^o^(O) > 
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which can be normalized by the known initial value of 
the coherence factor. Figure [T] presents the spectrum 
(^Sq){uj) and its Fourier transform {SQ{t)) for a wide 
range of external magnetic fields covering the full non- 
perturbative crossover Bfluc ^ B < A, where i?fluc = 
(Ej=i A'jy^'^ is a typical Overhauser field due to spin 
fluctuations [l8| , all the way down to very weak magnetic 
fields. All plots are obtained for an ensemble of = 36 
nuclear spins by sampling il = 10^ configurations. For 
any finite size system, the spectrums consist of a series 
of delta peaks which are smoothed into Lorentzians of 
width O.OOli? in the plots. We stress that this broadening 
is much smaller than the width of the peaks seen on the 
graphs, i.e. the peaks are composed of a large number of 
contributions at similar frequencies. 

The spectrum is basically characterized by two struc- 
tures. First we find a peak around the "bare" Larmor 
frequency B + given by the total effective magnetic 
field felt by the central spin. For the strongest mag- 
netic fields this sharp peak is the dominant feature whose 
nearly Lorentzian lineshape leads to exponentially decay- 
ing oscillations. As the magnetic field is lowered its width 
increases giving rise to faster decoherence. In addition, 
there is a low-frequency structure which carries a very 
low weight at strong fields but leads to slow modulations 
of the envelope function. These modulations were pre- 
viouslyobserved in the fourth order perturbative treat- 
ment [7| and arc now confirmed by our exact results. 

When lowering the field below Bfjuc ~ 0.1^ the low- 
frequency structure becomes the dominant feature, even- 
tually taking a scaling form of cj /B for sufficiently weak 
fields (see Fig. [T]c). The low- frequency structure has a 
finite width (oc B) and, as B — >■ 0, collapses into a delta 
peak. In the long-time dynamics, this leads to the slowly 
decaying low- frequency oscillations shown in Fig. [T]e), 
whose period and lifetime can be made arbitrarily large 
since the real-time evolution becomes a function of Bt. 
As i? — J> this ultimately gives rise to a non-decaying 
fraction representing nearly 0.5 of the initial value. At 
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FIG. 1. Top: Spectrum {S^{oj)) in the crossover and weak magnetic field regimes. The dashed lines mark the "bare" Larmor 
frequencies B + /ifnif Panel c) is plotted in terms of the rescaled frequency lu/B and shows only the low- frequency structure. 
Bottom: The corresponding real-time evolution of the coherence factor (^(^(t)). Each curve has an offset of 1 compared to the 
previous one. Black lines in panel d) are the norms [(^(^(t))! which, in the perturbative regime, correspond to the envelope 
function computed in 



B = the total coherence factor (Sq) + J2j'=i (^/) 
conserved, hence the lost central spin coherence is trans- 
ferred to the nuclear spins polarizing them along the cen- 
tral spin's initial orientation and ultimately locking the 
system into a non-decaying steady state. 

Here, by working with exact eigenstates of the model 
we get valuable insight in the processes involved which al- 
low us to demonstrate the complete absence of long-time 
decay, i.e. in contrast to other approaches [llj we do 
not observe a 1/ ln(t) decay. Indeed, we know that in the 
B ^ limit, solutions to the Bethe equations ([3]) split 
into two independent subsets. A state-dependent number 
< r < A/ of rapidities diverge as Xi « Li/B + 0(1) (Li 
being roots of a Laguerre polynomial For these, the 

creation operator S+(Ai oo) oc X]*;=o creates fully 
delocalized excitations whose contribution to the energy 
is 1/Ai cx _B -f 0{B'^). The remaining M — r rapidities 
stay finite Ai w A° -I- 0{B) and create localized excita- 
tions (due to the proximity of A,; with some particular 
values of the inverse couplings {e;}). They contribute a 
finite energy l/A,; oc const. -I- 0{B). These two subsets 
being decoupled, any set of finite rapidities solution 
to the B = Bethe equations ([3]) can be supplemented 
by any number r < N — rj of diverging rapidities to 
create exact B = eigenstates containing r + rf rapidi- 
ties, [S+(oo)]'' [UZi S+(A°)] ... which all have 
the same energy. This can be thought of as a manifesta- 



tion of BEC-like physics since it allows one to build eigen- 
states by adding any number of particles in a given single 
quantum state. It echoes the behavior of other integrable 
Gaudin models such as Dicke model's superradiance [2^ 
or Richardson model's superconductivity plj . 

Thus for any M-particle eigenstate |V'n), there exists 
a corresponding M + 1 eigenstate which, at B = 0, only 
differs by adding a single delocalized zero-energy parti- 
cle. Any such pair of eigenstates leads, in Eq. to 
a low-frequency contribution at w = ujm — Wn oc i? and 
contributes to the non-decaying fraction at i? = 0. Any 
other pair of eigenstates gives a contribution at a finite 
frequency which, being spread over a large energy band, 
leads to an initial rapid decay. Since the energies are only 
related to the finite rapidities content of the states, weak 
magnetic fields will only give subleading corrections. As 
shown in Fig. [2] for _B < -Bfiuc, this results in a satura- 
tion of the decay time characterizing the early dynamics 
which become well captured by a quadratic decay inde- 
pendent of the magnetic field value. 

Experimental observation of the non-perturbative 
long-time contributions, e.g. the low-frequency oscilla- 
tions, requires two conditions. It needs to have a lifetime 
sufficiently longer than the initial decay and, at the same 
time, it has to carry a sufficiently large weight. From 
a rough analysis of the real-time numerical data, fields 
B < OMA (~ 0.1 T for GaAs or ~ 0.0015T for Si:P 
using Q A/{g*HB) = 3.5 T and A/{g*HB) = 0.05 T 
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FIG. 2. Short time dynamics of the norm |(Sj(i))| for 
intermediate to weak fields. The nine values of magnetic 
fields B < -Bfluc shown in Fig. [T] are plotted with black lines 
{B/A e [1.3889 • 10"®, 0.05556]). The red line is a quadratic 
fit with parameter r ~ {A/N)~^. For both GaAs and Si:P, 
using values from 0] gives r ~ l^s. 



respectively) would be sufficient to exhibit oscillations 
with an amplitude of ^ 10% of the initial coherence fac- 
tor and a lifetime long enough compared to the initial 
rapid decay (see the B = 0.02778(^) curve in Fig. [T]e). 
These oscillations would be clearly distinguishable from 
any perturbative result since their frequency is nearly 
an order of magnitude lower than the Larmor frequency. 
Since we are looking for a clear measurable signature, 
the resulting threshold is at much weaker fields than the 
naive limit of validity of perturbation theory B A. 

In Fig. Owe show the long-time averaged contribution 
for a variety of system sizes. We see that even for mod- 
est system sizes {N > 20) the finite-size discretization 
does not seem to affect strongly the total weight car- 
ried by the low-frequency contributions. The error bars 
and system sizes treated here do not exclude a slow log- 
arithmic reduction ~ 1/ In {(3N) of the total fraction as 
was observed in a finite-size semi-classical treatment Q 
of the longitudinal decay, eventually leading to 1/ hi(t) 
decay in the continuum limit In contrast, our full 
quantum mechanical treatment includes the discreteness 
of the spin bath and shows the non-decaying fraction to 
be a observable for finite systems. Indeed, even an as- 
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FIG. 3. Non-decaying coherent fraction as a function of the 
number of nuclear spins at B/A = 1.04168 • 10~^. We sample 
Q — 10^ configurations. Error bars indicate the magnitude 
of the fluctuations which are due to a mixture of finite size 
effects and MC error. Dashed lines are guides to the eye at 
and 0.4812 (long-time values for A*' = 48). 



FIG. 4. Long-time averaged coherent fraction for 10 randomly 
generated initial nuclear spin configurations. The meaning of 
the error bars is as in Fig. O Inset: Corresponding real-time 
evolution. Magnetic field is B/A = 1.3889 x 10"®. 



sumed ~ 1/ In (/37V) reduction would yield an observable 
non-decaying fraction for experimentally relevant system 
sizes [N ^ 10^ in Si:P or TV « lO^^ in GaAs). 

Considering that even a physically narrowed system 
should be described in terms of a diagonal ensemble av- 
erage, in Fig. |4] we compare the weak field dynamics 
for an ensemble of 10 unpolarized {M ~ N/2) initial 
bath states obtained by choosing a random set of M up- 
pointing nuclear spins. This is less restrictive than the 
eigenstate content of a narrowed state since the initial 
Overhauser eigenvalue hf^^^. is not fixed. While we ob- 
serve slight variations in the total non-decaying fraction, 
it systematically stays above 0.4 and maintains the X-Y 
plane phase of the initial central spin orientation. Initial- 
izing the central spin at a different point on the Bloch 
sphere a\^) + /? I-U-) would only lead to a multiplicative 
factor {S^{t)) = a/3 (S^ {t)) which stiU conserves 
the relative phase. Albeit averaging may lead to a slightly 
lowered total fraction, the common phase obtained for 
all nuclear spin configurations means that diagonal en- 
semble averages (S^it)) = ^„ 7„ ({e}„| S^{t) |{e}„), be 
they performed for a narrowed or even infinite tempera- 
ture state, would still yield a large non-decaying fraction. 

Since it is linked to diverging rapidities, the existence 
of a " condensate" mode is completely independent of the 
set of couplings constants and hence the dot geometry. 
Therefore, the zero-frequency peak which leads to a non- 
decaying fraction, should not exhibit any inhomogeneous 
broadening in experiments which involve an average over 
different quantum dots. In fact, neither the geometry nor 
the initial nuclear spin configuration affect the resulting 
delta peak. This is in stark contrast to recent optical 
spin noise experiments (22j where a Lorentzian lineshape 
was observed. However, a direct comparison to our re- 
sults is not possible considering the experimental work 
was carried out on hole-spin based samples for which the 
hyperfine coupling is strongly anisotropic. 

In the isotropic system considered here a finite lifetime 
would require the addition of integrability breaking terms 
to the Hamiltonian ([T]). The weak dipolar coupling be- 
tween the nuclear spins could play such a role although 
on relatively long time scales Tdd « 10~* s in GaAs; an- 
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other mechanism would be quadrupole couphngs to the 
nuclear spins |23j . In addition, while we argued that en- 
semble fluctuations in the nuclear Overhauser field would 
not lead to a finite lifetime, it should be noted that local 
fluctuations in the external magnetic fleld would induce 
broadening due to the uj/B dependency of the spectrum. 
This demonstrates that, at weak field, the strong corre- 
lations between the central and nuclear spins lead to a 
clear distinction between external magnetic field and the 
internal Overhauser field. 

In conclusion, we have developed a non-perturbative 
method to study the dynamics in the central spin model 
which is applicable to arbitrary magnetic fields and ini- 
tial bath polarizations and provides access to the real- 
time dynamics at arbitrary long times. We applied this 
method to the transverse spin relaxation which, at weak 
fields, shows initial rapid quadratic decay followed by a 
large arbitrarily long-lived finite coherent fraction result- 
ing from the BEC-like behavior of the system. 
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